Synchronization of degree correlated physical networks 



in 
o 
o 

(N 
G 



C3 



I 

C 

O 
o 



> 
m 

O 

in 
o 

c3 



T3 

C 

O 
O 



x 



M. di Bernardo*, F. Garofalo*, F. Sorrentino*''' 

'Department of Systems and Computer Science 
University of Naples Federico II, Naples 80125, Italy 
' Corresponding author. E-mail: fsorrent@unina.it 
(February 2, 2008) 

We propose that negative degree correlation among nodes in a network of nonlinear oscillators, 
often detected in real world networks, is motivated by its positive effects on synchronizability. In 
so doing, we use a novel methodology to characterise degree correlation based on clustering the 
network vertices in p classes according to their degrees. Using this strategy we construct networks 
with an assigned power law distribution but changing degree correlation properties. We find that the 
network synchronizability improves when the network becomes dissassortative, i.e. when nodes with 
low degree are more likely to be connected to nodes with higher degree. Our numerical observations 
are confirmed by the analytical estimates found in this letter using an innovative approach based 
on the use of graph theoretic results. 



Networks of oscillators abound in physics, biology and 
engineering. Examples include communication networks, 
sensor networks, neuronal connectivity networks, biologic 
networks and food webs. Under certain conditions such 
networks are known to synchronize on a common evolu- 
tion with all the oscillators exhibiting the same asymp- 
totic trajectories. Moreover synchronization was ob- 
served to play an important role in many different prob- 
lems of a most diverse nature (physical, ecological, phys- 
iological, meteorological, to name a few); see for example 
[1-8]. 

Recently, it has been proposed that the network topol- 
ogy, i.e. the way in which the oscillators are mutually 
coupled between themselves, has an important effect on 
its synchronization properties. In [9], it was shown that 
the eigenratio R = between the highest eigenvalue 
An and the lowest eigenvalue A2 of the Laplacian associ- 
ated to the network structure is an essential measure of 
the network synchronizability, i.e. the smaller the eigen- 
ratio, the larger the interval of the values of the coupling 
gain, say a, for which the stability of the synchronous 
state is achieved. It is therefore important to charac- 
terise how the network topological features affects the 
Laplacian eigenratio. For example, scale free networks, 
which are common in nature, were found to show better 
synchronizability for increasing value of the power law 
exponent [10], [11]. 

Another important topological property of physical 
and biological networks is that often their nodes show 
preferential attachment to other nodes in the network 
according to their degree [12], [13]. According to this 
property, networks are said to exhibit assortative mixing 
(or positive correlation) if nodes of a given degree tend to 
be attached with higher likelihood to nodes with similar 
degree. (Similarly disassortative networks are those with 
nodes of higher degree more likely to be connected with 
nodes of lower degree.) 

The presence of correlation has been detected exper- 



imentally in many real-world networks. Interestingly, 
from the analysis of real networks, it was noticed that 
social networks are characterized by positive correlation, 
while physical and biological networks show typically a 
disassortative structure [12]. For instance, in [14], In- 
ternet was found to exhibit disassortative mixing at the 
AS level. A pressing open problem is to understand why 
negative correlation is an emerging property of physical 
and biological network. 

In this Letter, we study the effects of degree correlation 
on the network synchronizability properties both analyt- 
ically and numerically. Our main finding is that disassor- 
tative networks of oscillators synchronize better. These 
findings have immediate theoretical and experimental rel- 
evance for understanding the properties of real- world net- 
works. Because of its effects on synchronizability, we con- 
jecture that disassortative mixing has played the role of a 
self organizing principle in leading the formation of many 
physical networks as the Internet, the World Wide Web, 
protein interactions, neural and metabolic networks. 

We apply our analysis to scale-free networks in which 
the probability of finding a node having degree k scales as 
fc~ 7 . The network generation model is the configuration 
model as described in [15] [10]. The overall results are 
summarized in Fig. 1(a) where the effects of varying the 
degree correlation (measured through the index r defined 
below) on the Laplacian eigenratio R are shown for dif- 
ferent values of the degree distribution exponent 7. For 
all values of 7, we observe a reduction of R for decreas- 
ing values of r. This means that disassortative mixing 
enhances the network synchronizability. Interestingly, as 
depicted in Fig. 1(b) and Fig. 1(c), we observe that, un- 
der variations of the correlation parameter, the changes 
in R seem to be mainly due to variations of A2 while Xn, 
the largest eigenvalue of the Laplacian, is found to be 
practically independent from r. As discussed in [10], in 
the case of uncorrelated networks (r = 0), synchroniz- 
ability improves for increasing values of 7. 
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FIG. 1. Synchronizability of degree correlated scale free 
networks of size 10 3 nodes. Behavior of the eigenratio Aat/A2 
(a), of the second lowest eigenvalue A2 (b) and of the highest 
one Ajv (c), as functions of the correlation coefficient r defined 
in [11], for 7 varying from 2 (top line) to 5 (bottom line) in 
steps of 0.2. (d) The eigenratio as function of 7 as varying the 
correlation coefficient r from —0.3 (bottom line) to 0.3 (top 
line) in steps of 0.1. The lines are guided by the eye. 

We find that, as shown in Fig. 1(d), this is still the 
case when degree correlation is introduced. A consistent 
decrease of the values of R for the same value of 7 is 
observed when the network is disassortative, i.e. r < 0. 

Historically algebraic graph theory has been widely 
applied to describe structure and functions of networks. 
Here we will give evidence to the fact that tools from alge- 
braic graph theory can be conveniently used for the anal- 
ysis of network dynamical properties, such as the network 
synchronizability. The analytical framework we propose 
to capture the essence of the observed phenomenon is 
based on the following steps: (i) we introduce a novel 
index to estimate network assortativity and use such a 
quantity to construct networks with fixed degree distri- 
bution but different (desired) degree correlation proper- 
ties; (ii) using Cheeger inequalities from algebraic graph 
theory, we derive novel bounds on the Laplacian eigen- 
values of interest and hence new analytical bounds on the 
eigenratio R; (iii) using the new bounds, we validate our 
numerical observations. 

To start, we assume the network consists of N identi- 
cal oscillators coupled through the edges of the network 
itself [9], [10], [11]. Moreover, we suppose each oscillator 
is characterized by its own dynamics, {xi(t),i = 1...A}, 
described by a nonlinear vector field, say / = f(x), per- 
turbed by the output function of its neighbors repre- 
sented by another nonlinear term, say h = h(x). The 
equations of motion for each oscillator can then be given 
as follows: 



dxi 
~~dt 



N 



f(xi) -a^2£ijh(a 



1,2,. ..N, (1) 



where a is the overall coupling strength, and C = {Cij} 
is the Laplacian associated to the network topology [16]. 
Given a network of interest, we start by dividing the net- 
work vertices in p classes such that each one contains 



ni,n2—n p 



nodes and J2i 



N. In particular, we as- 



sign network vertices to each class in order of increasing 
(decreasing) degree. Hence, if we label as ki the mean de- 
gree of vertices belonging to the i-th class, then ki < fc,+i 
(i = l,2..,p-l). 

The probability, p(i), that a generic vertex belongs to 
class i is then given by p(i) = rii/N. Thus, considering 
all the vertices in each class i as they had mean degree ki, 
we can extend the usual definition of degree distribution, 
assuming that p(ki) = Ui/N . Analogously the probabil- 
ity of finding a vertex belonging to class i at the end of 
a randomly chosen edge within the network is given by 
q(i) = riiki/Y^iniki. 

Then, following [12], the presence of degree correla- 
tion with respect to the subdivision in p classes proposed 
above, can be estimated by using a new coefficient, r p 



defined as: 



k T (E qq T )k 



(2) 



where a q is the standard deviation of the distribution 
qk, k = ( ki, k 2 , k p ) T , q = ( qi, q 2 , q p ) T and E = 
{eij} € M pxp , with eij being the probability that a ran- 
domly chosen edge in the network connects nodes be- 
longing to class i and j (note that in the particular case 
where each class contains exactly only the vertices hav- 
ing a given degree, r p coincides with the coefficient r as 
defined in [12], [13].) 

From (2), it is possible to derive the distribution of 
edges among the network vertices as a function of r p as 
follows: 



E = qq J 



(3) 



where M is a symmetric matrix having all row sums 
equal to zero and appropriately normalized such that 
k T Mk = 1. Specifically, we can express M as follows: 



M 



mm 

(k T m) 2 ' 



where m = (mi, m 2 , m p ) is a vector such that 
X]i m i = (for instance we can choose m such that 
m, < rrii + i for i = 1, 2...,p — 1 in order to have a conve- 
nient form of the matrix M with positive values near the 
main diagonal and negative values far away from it). 

As explained in [17], it is possible to devise a strategy 
similar to the one presented in [12], [13], using the new 
coefficient r p , to generate networks with a given degree 
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distribution and a desired value of the degree correlation 
coefficient r p . Such a strategy was used to carry out the 
computations depicted in Fig. f . 

In [10] analytical bounds were given to explain the 
changes observed in the eigenratio R as the parameter 
7 was varied in a scale-free network topology. We found 
that, although these bounds should hold for any generic 
network topology, they seem to be inappropriate to ac- 
count for the changes in R observed in the network when 
degree correlation is introduced. In particular, as shown 
in Fig. 2, the values of the upper bounds on R computed 
according to the formulas given in [10] give estimates 
which do not reproduce the behavior of the eigenratio 
under variation of the network degree correlation. 

Therefore we shall seek to define new analytical bounds 
based on the mathematical theory of graph spectra. In 
particular, as Ajv was found to be almost independent 
from the correlation coefficient r (see Fig. 1(c)), we will 
focus on estimating the effects of correlation on the eigen- 
value A2, the parameter known as algebraic connectivity 
of graphs [16]. Specifically, given a graph, consider a sub- 
set of edges which disconnects the graph in two parts, also 
termed as a cut. For a given subset of vertices, say 5, we 
define hc(S) as the quantity given by: 

V(S)N 

hG{s) = \s\(N-\s\y (4) 

where T>(S) is the number of edges in the boundary of 5 
and |5| < , is the number of vertices in S. The Cheeger 
constant of a graph is given by ha = rains ho(S) and the 
so-called Cheeger inequality states that A2 < Hg [18]. 

We will show below that (4) can be successfully used 
to compute an upper bound on A2. To overcome the 
limitations due to the computation of the subset S that 
minimizes ha(S), we will follow a stochastic approach in 
order to estimate hc(S), starting from the available in- 
formation we have on the network. We will assume that 
the noticeable features of the network are only the de- 
gree distribution and the correlation specified; all other 
aspects being completely random. Note that finding the 
subset S such as to achieve the minimization of hc(S) is 
known to be an NP-hard problem [19]. 

Then, we can give a full characterization of a randomly 
chosen subset 5 in terms of the number of nodes in it, 
say Xi, belonging to each class i (i — 1,2..., p) and the 
network correlation measure r p . Analogously let us term 
as Vi — — the fraction of nodes in S drawn from each 
class i (i = l,2...,p). Note that the subset S is not 
supposed to satisfy any particular condition, not even of 
being connected. 

Now, we observe that the number of edges in the 
boundary, say T>(S), is given by the total number of edges 
starting from the vertices in 5, less the ones, say 1(5*), 
that are contained in 5, i.e. having both endpoints in 5. 
Thus we can estimate X and V as follows: 

l(S)=I(y,r p ) = £y T Ey, 

V{S) = V(y, r p ) = x T k - 21(5) = 2£ (y T q - y T Ey), 



where x = (xi,x 2 , ...,x p ) , y = (yi,2/2j — ,2/p) and £ 
is the total number of edges in the network. 
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FIG. 2. Upper and lower bounds given on the eigenratio 
Ajv/A2 as function of r (as defined in [11]) for different values 
of 7: 7 = 2.5(-) , 7 = 3(+), 7 = 4(0), 7 = 5(x), N = 10 3 . 
The lines are guided by the eye. The figure shows the compar- 
ison between the bounds in [2] (on the left) and those (9) (on 
the right). The behavior of the actual eigenratio is plotted in 
Fig. 1(a). 



Thus, hc{S) becomes: 

2£ (y 1 q - Y T Ey)N 



hc(y,r p ) 



(n T y)(N-n T y) ' 



(5) 



under the constraint that n T y < N/2, where n is the 
vector (ni,ri2, n p ) . A numerical optimization al- 
gorithm can then be used to find the subset S that 
minimizes /ic(y, r p ) in terms of y\, 2/2, Vp (and sub- 
sequently r p ) and, in turns, an upper bound for A 2 . 
Also, from (3) and (5), we get: 



dh G {S) k dViS) = _ 2£ )2 ^ Q 



dr„ 



dr n 



(6) 



Since, for any vector y, (6) is satisfied, then we have that 
< 0. Therefore, we can predict analytically that ha 
and hence A2 will be decreasing as the degree correlation 
is increased and, as a consequence, the degree eigenra- 
tio 4^ will increase for higher values of the correlation 
coefficient. 

As shown in Fig. 2, this is indeed what is observed 
with the new bound on A 2 giving a much better estimate 
of the behavior of the eigenratio with respect to both 
changes in the degree distribution and the degree corre- 
lation. (Note that A at is found to be almost independent 
from r in Fig. 1.) 
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FIG. 3. Effects of the maximization of A2 on the correlation 
parameter r defined in [11] for different scale-free networks 
with 7 = [2,3,4,5]. 

Another interesting inequality in spectral geometry is 
due to Mohar [19]: 



where h' G = min s ^ = 2£(y T q - y T £y)/(n T y). Us- 
ing (7), we can also get a lower bound on A2. Then fol- 
lowing an approach similar to the one used to compute 
the upper bound, it is easy to show that the lower bound 
in (7) has to decrease with r (note that when making the 
correlation change, the degree distribution is fixed and 
thus, k m ax cannot vary with r). Then since both the up- 
per and the lower bounds have to decrease with r, A2 is 
also expected to have the same trend. 

In order to derive bounds on R, we can use the novel 
upper and lower bounds on A2 and the bounds on Ajy 
(that was observed numerically not to depend on r) given 
in [20]: 

N 

~^k max < Ajv 5^ 2k max . (8) 

In so doing, we easily get the following analytical bounds 
on the Laplacian eigenratio: 

N k max Aat 2 ,qs 

N-l h G ~ A 2 " SSI" [ ) 

The comparison between these bounds computed as ex- 
plained above and those proposed in [10] is shown in Fig. 
2. We observe that the upper bounds in (9) provide bet- 
ter estimates of changes in the eigenratio under variations 
of 7 and more importantly r. 

The main result of the derivation presented above is 
the finding that disassortative networks synchronize bet- 
ter. Note that such a derivation can be easily extended to 
weighted topologies, i.e. to the case in which the strength 
of the coupling of each vertex on its neighbors is rescaled 



by the vertex degree. We wish now to assess whether 
negative correlation can be thought of as an emerging 
property of networks with an assigned degree distribution 
in order to improve their synchronization. Since Aat was 
found to be practically independent from variations of 
r, we use a simulated annealing meta-heuristic technique 
[21], to solve the problem of maximising A2 while keep- 
ing unchanged the network degree distribution. Namely, 
given a network with a certain degree distribution, we 
perform the following iterative procedure. At each itera- 
tion i, the endpoints of a randomly selected pair of edges 
are exchanged if exp ( ~ A j, A2 - ) ) > z; z being an uniformly 
distributed random variable between and 1, A(A 2 ) the 
variation achieved in the objective function A2 before and 
after the execution of the move and T the system tem- 
perature. 

As shown in Fig. 3, under the effects of this iterative 
procedure aimed at maximising A2 , we observe the spon- 
taneous emergence in the network of interest of negative 
degree correlation, namely the increasing of its disassor- 
tativity. Thus, following an entirely different approach, 
we come to the conclusion that if the degree distribution 
of a given network is fixed, then to improve its synchro- 
nizability one has to introduce negative degree correla- 
tion among its nodes. 

In conclusion, we have shown that degree correlation 
among the nodes of a network of nonlinear oscillators 
does indeed have an effect on the Laplacian eigenratio 
and hence its synchronizability. Using a novel correlation 
index, we were able to show both analytically and numer- 
ically that the eigenratio is lower in disassortative net- 
works. Hence we conjectured that in many evolutionary 
physical and biological networks of nonlinear oscillators, 
negative correlation among nodes can be an emerging 
property aimed at facilitating the synchronization pro- 
cess. An iterative nonlinear optimization technique was 
used to illustrate this latter point. 
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